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A lattice gas model for Schlogl's second chemical reaction is described and analyzed. 
, Because the lattice gas does not obey a semi-detailed-balance condition, the equilibria are 

f***) ' non-Gibbsian. In spite of this, a self-consistent set of equations for the exact homogeneous 

equilibria are described, using a generalized cluster-expansion scheme. These equations are 
solved in the two-particle BBGKY approximation, and the results are compared to numerical 
' experiment. It is found that this approximation describes the equilibria far more accurately 

ON , than the Boltzmann approximation. It is also found, however, that spurious solutions to 

""^5 ' the equilibrium equations appear which can only be removed by including effects due to 

' three-particle correlations. 

wo: 
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I. INTRODUCTION 

Lattice gas automata have been widely used as models of nonequililprium statistical systems since it 
was shown in 1986 that they could be used to model Navier-Stokes fiuidsEl. Lattice gases consist of particles 
moving about and colliding on a lattice in such a way that their macroscopic behavior satisfies hydrodynamic 
partial differential equations. Like the Ising model, they are simple discrete systems which are well suited 
both to computer implementation and to elegant analytic techniques; unlike the Ising model, however, they 
can be used to study phenomena far from equilibrium. 

All of the usual tools of kinetic theory can be used for the analysis of lattice gases. Lattice gases whose 
collisions obey a condition known as semi- detailed balance (SDB) can be showno to have a Gibbsian (product) 
equilibrium distribution. As the lattice spacing goes to zero, expansion about this equilibrium yields the 
hydrodynapiic equations satisfied by the system; this is a discrete version of the usual Chapman-Enskog 
procedure^. 



*This work was supported in part by the divisions of Applied Mathematics of the U.S. Department of Energy 
(DOE) under contracts DE-FG02-88ER25065 and DE-FG02-88ER25066, and in part by the U.S. Department of 
Energy (DOE) under cooperative agreement DE-FC02-94ER40818. 
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To date, most analyses of lattice gases have been done using the Boltzmann molecular chaos approximation. 
We have recently used cluster expansion methods to develop an exact description of SDB lattice gasesu. In 
such lattice gases, the exact equations of motion include the effects of correlations which renormalize the 
lattice gas transport coefficients. In this paper, we extend these methods to describe a particular non-semi- 
dctailed-balance (NSDB) lattice rgas. Related work on exact equations for NSDB lattice gases has recently 
been done by Bussemaker et. al.B. 

It has been known for decades that chemically reacting systems far from equilibrium can exhibit fascinating 
phenomenology, including pattern formations and symmetry breaking!! Such complicated phenomenology 
can arise from very simple chemical reactions, and idealized model reactions have been developed to illustrate 
these phenomena. For example, the simple model reaction proposed by Schlogl in 1972cl, 

2X + A r± 3X, 

where X is the reactant species and A is a background species of fixed density, can posess two stable 
equilibrium concentrations of the species X. In that case, the system can exhibit spontaneous pattern 
formation as it breaks into domains of each concentration. Because kinetic fluctuations are important in 
the dynamics of such systems, it is natural that lattica-ga s..a utomata be applied to their study, and this has 
been done with great success over the past five years&EIOij. 

Reaction-diffusion lattice gas models typically allow reactant particles to diffuse for some number of 
timesteps, k, between reactions. The diffusion steps obey SDB, while the reaction steps usually do not. 
It is remarkable that while natural chemically reacting systems seem to be able to spontaneously generate 
patterns with microscopically reversible laws of motion, all lattice gas models of such systems to date have 
found it necessary to violate SDB. There is no doubt that it is easier to generate nontrivial structure in NSDB 
lattice gases. Violations of SDB can lead to the spontaneous generation of patterns and correlations, and 
hence non-Gibbsian equilibriatJ. In such situations, however, the Boltzmann molecular chaos assumption 
is particularly suspect, and the theoretical analysis of the system becomes difficult or impossible. Only in 
the limit of largeri; has analytic progress been made to date; at low k the Boltzmann theory is known to be 
seriously in erroiO. 

In this paper, we describe a simple lattice gas model for Schlogl's second chemical reaction. Because the 
reaction steps of this lattice gas do not obey SDB, the equilibria are non-Gibbsian. We derive a self-consistent 
set of equations for the exact homogeneous equilibria using cluster-expansion methods. We solve these 
equations in the two-particle BBGKY approximation; in this approximation these equations are similar to 
those arising from the method recently developed by Bussemaker et. aljj. Comparing our results to numerical 
experiment, we find that this approximation describes the equilibria far more accurately than the Boltzmann 
approximation. We also find, however, that spurious solutions to the equilibrium equations appear which 
can only be removed by including effects due to three-particle correlations. These spurious solutions are an 
important artifact of this technique, and we argue that it is necessary to pay very close attention to them in 
any such analysis. 



II. DESCRIPTION OF THE SCHLOGL MODEL LATTICE GAS 



A. Schlogl's Second Chemical Reaction 

Our starting point is the following generalization of Schlogl's second chemical reaction!: 

2X + A % 3X 
X + B^2X 
C^X, 

where X is the reactant species, A, B, and C are background species of fixed density, and the k^ are the 
forward (+) and reverse (— ) rates for the reaction with j reactant molecules on the left. Denoting the density 
of species Y by Ny, the stoichiometric equation for this reaction is 
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'' X X k+N A N x - k 2 N x + k+N B N x - KN 2 X + k+N c - k^N x 
k - k x N x + n 2 N x - n 3 N x 



dt 



where we have defined the stochiometric coefficients, 

k = k^N c 

«i = &o — k^N B 

k 2 = k^NA- fcf 

K 3 = fe J . 

Finally, to model the stochastic motion of the reactant X between reactions, we add a diffusive term to 
obtain the reaction-diffusion equation, 

^ =D\7 2 N x + Ko - Kl N x + k 2 N x -k 3 N x . (f) 

Note that Eq. ([!]) allows for up to three spatially uniform equilibria, corresponding to the roots of the 
cubic. When there are three roots and k 3 > 0, the low-density and high-density roots, denoted by N x and 
N x respectively, are easily seen to be stable to small fluctuations, while the middle root, N x , is unstable. 
The evolution of Eq. (|l|) from generic initial conditions thus yields domains of constant density N x and N x , 
separated by sharp gradients whose widths are governed by the diffusive term in Eq. (Q). (See Fig. ||.) 



B. Lattice Gas Model 



We model the kinetics of the generalized Schlogl reaction by a lattice gas automaton. This consists of a 
regular lattice, C, with n lattice vectors at each site; we denote the lattice vectors by c*, where i € {1, . . . , n). 
The state of the system at time t is then completely specified by the quantities n*(x, t) s {0,f} where 
i € {1, . . . , n} and x 6 C. We have n l (x, t) = 1 if there is a particle with velocity Ci at position x at time t, 
and n l (x, t) = otherwise. 

The evolution of the lattice gas for one timestep takes place in two substeps. In the propagation substep, 
the particles simply move along their corresponding lattice vectors, 

n*(x + Ci,t + At) <- n 4 (x, t). 

This is followed by the collision substep, in which the newly arrived particles change their state. The 
collisions are chosen to model the reactive and diffusive dynamics of species X. Their effect is captured in 
the collision operator, lo 1 , which gives the increase in the number of particles moving along direction i due to 
collisions. In terms of this collision operator, the full equation of evolution of the lattice gas may be written 

n l (x + c l ,t + At) = n\x,t) + uj 1 (n* (x, t)) , (2) 

where the dependence of uj 1 on n* (x, t) indicates that each component of the collision operator can depend 
on all the components n l at the local site. 

In this work, we restrict our attention to the Schlogl model in two dimensions. We use a hexagonal 
(honeycomb) lattice because it has only three bits of state at each site (n — 3), thereby greatly simplifying 
the analysis; at the same time, it is sufficiently symmetric to ensure the isotropic form of the density balance 
equation, Eq. (|f|). This lattice is illustrated in Fig. ([!]). Note that such a lattice can be colored like a 
checkerboard; note also that the correspondence between the bits and the lattice vectors is rotated by tt/3 
for the differently colored sites. 
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C. The Collision Operator 



Following previous work on the modeling of chemical reactions by lattice gascsElBEle, we define two 
types of interparticle collisions. The chemical reactions take place in reactive collisions in which particle 
number does not need to be conserved. Between reactions, the particles execute diffusive collisions in which 
particle number is conserved. Both types of collision processes are stochastic; that is, the outgoing state of a 
collision depends on one or more random bits that must be generated at each site at each time step, as well 
as on the incoming state. Reactive collisions occur once every k timestcps; the remainder of the collisions 
are diffusive. 

We need to carefully define the dynamics of the reactive and diffusive collisions, and thence the form of 
the respective collision operators, u>r and cJq. Because there are three bits per site, each site can be in 
one of eight states. We enumerate these states by specifying the three bit values, i.e., 000, 001, . . . 111. The 
collision process can then be completely determined by specifying the outgoing state corresponding to each 
incoming state. Since the lattice gas is stochastic, this specification may depend on one or more random 
bits. 

Let a(s — > s') be 1 if a collision takes state s to state s', and otherwise. Clearly, for each incoming state 
s, a(s —>■ s') can equal 1 for exactly one s', and must equal for all others. In terms of this transition matrix, 
the collision operator can be written 

n 

wV) = -» s ')( sH - s *) n <w > (3) 

s,s' j=l 

where = 1 — i — j + 2ij is the Kronecker delta of the two bits i and j. Together, Eqs. (||) and (|J) are a 
complete specification of the dynamics of the lattice gas. Note that a(s — ► s') may depend on random bits. 

D. The Boltzmann Equation 

We now suppose that we have prepared an ensemble of lattice gas simulations, on grids of the same size, 
with initial conditions that are sampled from some distribution. We may take averages across this ensemble. 
Denoting N l (x, t) = (n ! (x, t)}, the ensemble average of Eq. (^) is 

N l {x + Cl ,t + At) = N l (x,t) + (w*(n*(x,t))}. 

We are hampered from taking the ensemble average of the collision operator, Eq. (^), by the fact that it 
is generally a nonlinear function of the n l (x., t), and the average of the product is not equal to the product 
of the averages unless the quantities involved are uncorrelated. The simplest approximation to make is the 
Boltzmann molecular chaos assumption that the particles entering a collision are uncorrelated; in this case, 
the ensemble average of uj 1 yields the Boltzmann collision operator, 

n 

s,s' j=l 

where A(s — » s') = (a(s — ► s')) G [0, 1] is the ensemble-averaged transition matrix. 

Note that there are three one-particle states (001, 010, 100), three two-particle states (110, 101, 011), one 
zero-particle state (000), and one three-particle state (111). Let \s\ denote the number of particles in state s, 
so for example 1 1 1 1 = 2. For the lattice gas considered here, the mean outcome of both diffusive and reactive 
collisions depends only on the total number of incoming particles, and is always uniformly distributed over 
the states of the outgoing particle number. Mathematically, this means that the A(s — > s') can depend only 
on \s\ and |s'|, and can thus be tabulated as in Table |, where P*j is the probability that a collision will take 
a state with j particles into a state with i particles. The evolution equation in this approximation is the 
Boltzmann equation, 
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AT i (x + c i ,i + At) = N^^t) + O i (n*(x,£)). (4) 
For the diffusive collisions, we must have 

p\ = §*. 

3 3 5 

where J 1 ■ is the Kronecker delta. Thus, a diffusive collision is nothing more than a random permutation of 
the three incoming bits. Calculation of the corresponding Boltzmann collision operator is straightforward 
yielding 



Vt l D {N* 



-N l 



1 



-N 



i+l 



1 



-N 



i+2 



(5) 



where the superscript of N is understood to be taken modulo 3. 

To simplify the algebra for the reaction step, we henceforth restrict our attention to the following specific 
values for the particle transition probabilities, 



P' 



/2/3 1/3 

2/3 1/3 

1/3 2/3 

\ 1/3 2/3 



for i,j £ {0, 1, 2, 3}. Calculation of the corresponding Boltzmann collision operator yields 

f4(iv*) 



- - N l + - (^N 1 + N°N 2 + N X N 2 ) - —^N 1 ^. 
9 9 v ' 9 



(6) 



A complete Boltzmann description of the system is given by Eq. (||), using Eq. (g|) once every k timesteps 
and Eq. (||) otherwise. 



E. Boltzmann Equilibria 

Note that the Boltzmann equation, Eq. ([|), admits homogeneous, isotropic equilibria, iV = iV 1 = TV 2 = /, 
where / obeys f2(/) = 0. Note also that the diffusive collision operator, Eq. satisfies Q,jj{f) = 
identically. We thus find homogeneous, isotropic equilibria by demanding that the reaction step do likewise, 

- 

9 3 9 J 

= i(l-2/)(7/ 2 -7/ + l). (7) 

This has roots at / = | and / = \ {\ ± yf}' ® displays the evolution of the lattice gas model for 
these parameters, with the initial condition f = \ everywhere. 



III. EXACT EQUATIONS OF MOTION 



The exact microscopic equations of motion for any lattice gas-.are easily described in terms of the multi- 
particle means N a (following the notation of our previous papeicl we denote by a an arbitrary subset of the 
bits (particles) in the system, and by N a the ensemble average of the product of those bits). In terms of 
these means, the exact time-development equation is 
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N a (t + At) = A a K^(t), (8) 

where we use the convention of summing over any index which appears twice on one side of an equation and 
not at all on the other side. In this equation, A a ^ is an advection operator, described by a permutation 
matrix on the set of bit sets a, which carries each bit of the system forward along its associated velocity 
vector. The operator describes the collision process. It can be factorized into contributions from each 
lattice site, 

K\ = J] V\, (9) 

where Lp is the set of vertices associated with bits in (3 and /3 X is the set of bits in /3 at the lattice site x. 
The mean vertex coefficients V are related to the state transition probabilities A(s — > s') through 

n = EEH) HHsU (^ s ')' ( 10 ) 

The exact time-development equation (^|) can be rewritten in terms of connected correlation functions 
(CCF's) using the standard cluster expansion. The means are expressed in terms of the CCF's through 

N a = f a (T*)= r Ci r f2 ...r c ', (11) 

(eTr(a) 

where n(a) is the set of all partitions of a into disjoint subsets, Ci, ■ ■ ■ > Cq- F° r example, we have N a = T a , 
N ab = T ab + r apb Thig relation can be i nverte d to express the CCF's in terms of the means, T a = g a (N*). 
We can now rewrite (||) as 

T a (t + At) = A a g (K\P(r*)), (12) 

This exact equation has been used as a starting point in previous worksSS. It has been applied to SDB 
lattice gasesu, where the equilibria have no correlations and the expression on the right hand side can be 
linearized in terms of, the CCF's T a with \a\ > 2. Eq. ( |l2|) has also been applied to NSDB lattice gases 
by Bussemaker et al.fl who neglected CCF's T a with \a\ > 3, and thereby derived the 2-particle BBGKY 
equations for NSDB lattice gases. 

It has been showncl that the linearized form of ( |l2| ) can naturally be expressed in terms of a sum over 
diagrams, each of which is weighted by a product of factors associated with each vertex at each time step. 
There arc a finite number of possible vertices, so that a complete formulation of the dynamics of a SDB 
lattice gas can be given in terms of a set of "Feynman rules" for allowed diagrams and vertex weights. 

An analogous diagrammatic description can be given for the exact nonlinear equations ( |l2| ) . The nonlinear 
diagrammatic expansion can be derived by proving a general factorization theorem for the time development 
of CCF's including particles at different vertices. The essential ingredient in proving this factorization is the 
observation that if a set of variables a depends stochastically on another set of variables (3, so that the CCF 
T Q is given by 

T" = )C a 6 Y[ T*' (13) 

where £ = {£i, . . . ,£ m } is a set of (not necessarily disjoint) subsets of /3, then the CCF of a joined with a 
set of variables 7 which are not dependent on (3 is given by 

r u i r = jc« £ jjr«* uc * (14) 

C6^ m (7) « 

where £ = {(1, ■ ■ ■ , Cm} is summed over all partitions of 7 into precisely m distinct sets. This result essentially 
states that once we know an expression for the outgoing CCF's at a particular vertex of a lattice gas in terms 
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of the incoming CCF's, we can calculate the CCF of a set of particles at multiple lattice sites by applying 
( |l4| ) at each vertex separately. The general expression for an outgoing CCF at one vertex can be written by 
expanding 

$"(r*) = g^K^PiT*)). 
as an explicit polynomial in the CCF's; i.e., 

*'(n=K* n 1 * ( 15 ) 

where £ = {£1, . . . is summed over all sets of CCF's with nonzero coefficients. Each time the equation 
( |l4| ) is applied at a particular vertex, the correlated quantities at the other vertices are carried along and 
divided up in all possible ways among the incoming CCF's. A simple example of this result is that when a 
is an outgoing particle from a vertex with incoming particles 61,62,^3, and c is an outgoing particle from a 
different vertex at the same time step, we have (for a general lattice gas) 

and 

r ac _ &f r bjc 1 I — r b ' biC -+- — r bib2&3C n ri 

i ijij 

The proof of (Q) follows fairly easily by induction on j and k. The details of this proof and the general 
factorization theorem in the nonlinear case will be given in a separate publicationij. The result (|l6|), which 
follows directly from (|ll]) will be sufficient for our purposes in this paper. 

We conclude this section with a derivation of a simple form of the factorization theorem which we will 
need in the sequel. Assume that at one vertex we have an outgoing particle A and incoming particles a, 6, c, 
and that at another vertex we have an outgoing particle A and incoming particles a, b, c. We wish to find the 
dependence of the outgoing CCF T AA on the incoming correlations, neglecting all CCF's between 3 or more 
particles. It will suffice for us to know the dependence of the outgoing 1-particle means on the incoming 1- 
and 2- particle CCF's at each vertex. Thus, we can write 

T A = f( T a^ p 0j pc^ pafc ; pbc ; pac) + 

and 

r A = g{r a , r b , r e , r sb , r bs , r se ) + o(c a ) 

where by 0(C,) we denote quantities dependent on CCF's of i or more variables. Applying (|l6| ) once, we 
have 



r -&r« T + af^ r + 9f^ r +0{C3) - 



Applying ( |T^ ) again, we have 



^ AA = §^-^E^ + 0(C 3 ), (17) 



where a, a are summed over {a,b,c} and {a, 6, c} respectively. Note that this equation has a diagrammatic 
interpretation because the coefficient associated with the propagation of a pair of correlated quantities at 
different vertices factorizes into contributions from each vertex separately. We will use this simple factoriza- 
tion result in the next section to compute the exact 2-particlc BBGKY equations for the equilibria of the 
Schlogl model lattice gas. 
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IV. EXACT EQUILIBRIA OF SCHLOGL MODEL 



We will now consider the exact equations of motion for the Schlogl model lattice gas defined in Section 2. 
By neglecting correlations between more than two particles, we arrive at the 2-particle BBGKY equations, 
which we then solve using the diagrammatic method. The 2-particlc BBGKY equations were described 
for a general NSDB lattice gas by Bussemaker et al.El, who gave an iterative method for finding solutions 
to these equations. Although the equations we are solving here are essentially equivalent to those which 
would be found by applying the methods of these authors to the Schlogl model lattice gas, our diagrammatic 
method of solution of these equations is rather different. Using the diagrammatic formalism, there is no 
issue of convergence as there is with the iterative method; furthermore, in our analysis, there is no question 
of uniqueness of solutions - we can identify directly all distinct solutions of the 2-particle equations. In fact, 
we find that the 2-particle BBGKY equations have spurious solutions for the lattice gas considered here. 

The first step in writing the exact equations for the Schlogl model lattice gas is to write the exact equation 
for CCF's at a single vertex. There are two sets of such equations, corresponding to the diffusive and reactive 
vertices, respectively. The mean vertex coefficients V a p for both of these vertex types are symmetric with 
respect to permutations of incoming and outgoing bits separately, and therefore are only functions of the 
numbers of bits in a and (3. These vertex coefficients are easily calculated and are tabulated in Tables |H| 
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From these vertex coefficients, we can use (|12|) to write the exact equations for the outgoing CCF's from 
a diffusive or reactive vertex in terms of the incoming CCF's. These equations are again invariant under 
arbitrary independent permutations of the incoming and outgoing bits. Labeling the outgoing particles by 
A, B, C and the incoming particles by a, 6, c, the equations for a diffusive vertex are given by 

T A _ +r b + r c) 

3 

r AB = i(r ab + r ac + r bc ) (18) 

o 

The 1-particle equation for a reactive vertex is 

17 

pA _|_ ^papb _|_ papc _|_ pbpc _ 2p a p^p c _|_ pib 2p c T^ a ^ (19) 

9 9 

+r a C _ 2 r b r ac + r bc - 2r a r bc - 2r abc ). (20) 

The equations for 2- and 3-particle outgoing CCF's are straightforward to calculate but are algebraically more 
complicated than Eq. (^0|) . Note that setting the two- and three-particle correlations to zero in this equation, 
and setting all 1-particle correlations to the mean occupation number / = T a = T b = T c , reproduces the 
Boltzmann equilibrium, Eq. (j^). 

Henceforth, we will restrict attention to uniform equilibria, so that the correlations are independent of 
spatial coordinate or orientation. We denote the equilibrium values of the 1-, 2-, and 3-particle CCF's 
entering a reactive vertex by I\, I2, and -Z3 respectively. Similarly, we denote the CCF's leaving a reactive 
vertex by 0\, O2, and O3. The exact equations of motion for the 1- and 2-particle CCF's leaving a reactive 
vertex are 

(h- - — 



Q — I 1 _ 1 _ 1 4. ^ 1 _ tlZll 4. 2 _ _ -^l-^ (21) 

2 111 ~ ri-7 Q1 n 0-7 Q1 1-7 1-7 fl \ ' 



1 71? 


14/3 7h UIih 14Ja 




9 + X- 


9 3 3 9 




-1 


98/j 5 49/j 1 196/f 196/f 49/ 2 


98/i/ 2 98/^/2 


81 + 27 


81 9 ' 27 81 ' 27 


27 9 


784/3/2 


392/4/2 49/| 196/i/f 196/f/| 


98/ 3 196/f/ 3 392/3/3 


27 


27 9 ' 9 9 


81 ' 27 81 


196/ 2 / 3 


392/1/2/3 196/f 




27 


27 81 
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The equation for O3 can be similarly written, but is slightly more complicated and will not be used here. 
Recall that, as was demonstrated in the previous section, the exact dynamical equation of an arbitrary 
number of correlated quantities can be described in terms of the exact equations for the CCF's at a single 
vertex. Thus, Eq. (|2l]), along with the corresponding equation for O3, gives a complete description of the 
equations of motion of all CCF's at a reactive timestep. 



To complete the equilibrium equations (21), we must determine the relations between the outgoing corre- 
lations Oi from a reactive vertex and the incoming correlations Ii. Referring back to Eqs. ([l2]) and (]l8|), we 
see that at diffusive time steps, the correlations essentially perform random walks on the honeycomb lattice. 
Thus, the correlation I± entering a fixed reactive vertex at some time step is a weighted sum of outgoing 
correlations 0\ from vertices at the previous reactive timestep, with total weight 1. Since we have assumed 
an isotropic equilibrium, we have an equilibrium density / satisfying 

f = h=O v (22) 

It is interesting to note that by using this equality in the first equation in ( ^l| ) we can write an exact 
expression for ijj in terms of / 2 and I\. Inserting this expression into the second equation of (pl|), we find 
that the terms in / 2 cancel and we have the result 

02 = -\ + f-f (23) 

Note that this equation is exact, and must be satisfied by any isotropic equilibrium of the system. 

In principle, we would now like to find an exact set of expressions relating the quantities I2 , I3 to outgoing 
quantities 2 , O3 by iterating the exact equations of motion. However, this is technically infeasible since such 
a calculation would involve a sum over diagrams involving arbitrary numbers of correlated quantities. Thus, 
we shall now restrict to the 2-particle BBGKY equations by neglecting correlations of more than 2 particles. 
By making this simplification, we derive a simple set of equations whose solutions give the equilibria of the 
lattice gas in the 2-particle BBGKY approximation. 

Neglecting 3-particle correlations, and setting f = 0± = Ii, the exact equations for the CCF's at a reactive 
vertex become 

1 , 7/ 2 14/ 3 , 1I 2 14// 2 



/ = 



9 3 9 3 3 
1 2 



= / + -(!- 2/)(l - 7/ + 7f z + 21 J 2 ) (24) 



-1 49/ 2 98/ 3 49/ 4 196/ 5 196/ 6 49/ 2 

O2 



(25) 



81 27 81 9 27 81 27 

98/J 2 98/ 2 / 2 784/ 3 / 2 392/ 4 / 2 49/ 2 2 196/Jf 196/ 2 / 2 
27 9 + 27 27 9~~ + I) 9 ' 

The first of these equations is satisfied whenever either 
or 

^ 2 = -^(l-7/ + 7/ 2 ) (26) 

The solution / = 1/2 corresponds to the unstable equilibrium of the Boltzmann theory, and shows that this 
unstable equilibrium still exists in the 2-particle BBGKY approximation. We will not discuss this solution 
further here. Inserting (|2^) into (p5|), we again derive the identity (|3|), so this identity still holds in the 
2-particle BBGKY approximation. 
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To find all solutions to the 2-particle BBGKY equilibrium equations, it remains for us to find a relation 
between I2 and O2 ■ The analysis of the flow of 2-particle correlations is slightly more subtle than that of the 
1-particle density. Tracing back a given incoming correlation I2 to the previous reactive vertex (k timesteps 
earlier), we find that with some probability ^fc(l) the random walks of the correlated quantities lead back 
to a pair of outgoing particles from a single vertex associated with an outgoing correlation 2 - However, 
the remaining random walks (with probability 1 — </>/c(l)) lead to a pair of correlated quantities at different 
vertices. For a fixed pair of vertices, we denote such an outgoing correlation from a reactive vertex by 0\ t \. 
Using Eqs. ( pp[ ) and ([l7|), we can expand 0\.\ in terms of incoming CCF's of the 6 particles associated with 
the two vertices in question. Making the 2-particle BBGKY approximation, we have 

M =A 2 (£l M ) (27) 

pairs 

where the sum is taken over all 9 possible pairs of incoming particles, 1 from each vertex, which may be 
correlated, and where 

A = A(/,7 2 ) = y (/ - f - I 2 ) = 1(1 + 14/ - 14/ 2 ). (28) 

Note that A < 1/3, with equality only when / = 1/2. We can repeat the above steps for the particles 
correlated in each term Moving back through k — 1 diffusive vertices, associated with random walks of 
the correlated quantities, we again have some set of diagrams where the correlation originates in a pair of 
outgoing particles from a single previous reactive vertex, and some other set of diagrams where the correlated 
quantities are still separate. Repeating this analysis indefinitely, we find that the equilibrium correlations I2 
and O2 can be related by 

00 

/ 2 = []T<Mi)(3A) 2t - 2 ]0 2 , (29) 
t=i 

where <j>k(t) is the weighted sum over all diagrams describing random walks of 2 particles for kt time steps 
on the honeycomb lattice, where the particles leave a particular vertex on the first step in a fixed pair of 
directions and arrive together at some possibly different vertex at the final step. In these diagrams, the 
particles are not allowed to visit the same vertex at any time step divisible by k (reactive vertices), and when 
they visit the same vertex at any other time step (diffusive vertices), they exit in different directions with 
each possible pair of outgoing directions having equal probability (corresponding to (|l8|)). Note that the 
factor of 3 appears because the usual probability 1/3 of a given random bounce is replaced by the weight A. 
As an example of a coefficient <f>k (t) , it is easy to calculate 

Mi) = 1, 

since the unique diagram which contributes is as shown in Figure [?]. Similarly, since at t = 2 there are 30 
diagrams which each contribute (1/3) 6 , one finds that 

<■»-» (*)'-£■ 

It follows immediately from the random walk interpretation of (j>k(t) that 

00 

X>(*) = x > 
t=i 

since the probability that two random walkers in 2D will eventually collide is 1 . An immediate consequence 
is that the series 
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X>(*)(3A) S 



t=i 



converges whenever A < 1/3. Furthermore, for A satisfying this condition, we can calculate the above scries 
to arbitrary accuracy; gi 1 
follows immediately that 



to arbitrary accuracy; given any e, when A < 1/3 we can choose T such that 4>k{t) > 1 — e, and it 



5>*(*)(3A) 



2t - 2 < e. 



Thus, to calculate the sum to within an accuracy of e we need only calculate a finite number of coefficients 
<pk{t), a task which is easily performed numerically by a computer. 

We may now use (|23|), (p6[), and ( p9| ) to derive a single equation for the 2-particle BBGKY equilibrium 
density /, 

C(/) = 3(1 - 7/ + If) - 7(1 - 9/ + 9/>(/) = 0, (30) 

where 

"2 



«(/) = !>*(*) 



*=i 



g (l + 14/-14/ 2 ) 



(31) 



Since, as mentioned above, we can calculate a(/) to an arbitrary degree of accuracy, it is a straightforward 
process to numerically determine the values of / which satisfy Eq. (^50|) to an arbitrary degree of accuracy. 
We have performed such a numerical analysis for k ranging from 2 to 7. For each value of fc, we find not 
two, but four distinct equilibria satisfying C(/) = 0- As an example, we graph in Figure ^ the function £(/) 
for k — 3. This function has two zeros at / ~ 0.1903 and / 0.8097, which presumably correspond to the 
actual equilibria of the system. We will refer to these zeros as the "primary" solutions. In addition, however, 
the function has two zeros near / = 1/2, which we will call the "secondary" solutions. Because the series for 
a(f) converges very slowly in the vicinity of / = 1/2, one might be suspicious of the secondary solutions. 
To see that such solutions must exist, however, we can observe that at / = 1/2 we have a(l/2) = 1 and 
therefore C(l/2) = 13/2 for any k. Since the function converges nicely and is negative above (below) the 
lower (upper) primary solution, there must be a secondary pair of solutions, just as we see in the graph. 

Comparison of the primary equilibrium solutions with numerical results from simulations of the lattice 
gas with various values of k shows that these solutions of the 2-particle BBGKY equations predict the exact 
equilibria of the lattice gas system remarkably well. This comparison is given in Fig. ||. We see that the 
2-particle BBGKY approximation gives an excellent numerical prediction of the equilibria of the Schlogl 
model lattice gas. However, the existence of the spurious secondary equilibria demonstrates emphatically 
that one must be very careful when dealing with truncations of the exact equations for a lattice gas. In 
the workppf Bussemaker et al., for instance, an iterative method is used to solve the 2-particle BBGKY 
equationsa. This approach can result in a spurious equilibrium, with no indication that any other solution 
exists. Thus, without some further criterion for judging the validity of a solution to these equations, it is 
difficult to evaluate the results of such an analysis. 

We will now proceed to give some simple analytic arguments which show that the secondary solutions 
are highly sensitive to the introduction of 3-particle CCF's, and thus that they are suspect from a priori 
grounds. First, let us observe that the introduction of a small amount of 3-particle correlation in I3 would 
change (^6|), which would then read 

* -4(1 -*/ + */"> + 3(1^7)- < 32 > 

If the correlation I3 were small, this would cause a change in I2 which would be small except in the region 
/ ~ 1/2, where the change would be dramatic. A change in I2 would in turn cause a comparable change in 
A through (p8[). Since the sum ([n]) converges slowly in the region of A « 1/3, the value of a(f) is highly 
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sensitive to a slight change in A in this region, which is precisely the region where / « 1/2. In fact, only a 
small change in A is needed to lower a sufficiently that £(1/2) < 0, which would result in a disappearance of 
the spurious equilibria. 

The composition of the two extreme sensitivities described here makes it clear that the existence of the 
spurious equilibria arc highly dependent upon the vanishing of the 3-particle CCF I3. In fact, we have 
extended our analysis to include a simple class of 3-particle diagrams and found that with this minor 
modification, the spurious equilibria completely disappear. Specifically, one can take the exact 3-particle 
equations at a vertex, and solve using the additional condition that I3 = /XO3 where \i is the weight of some 
simple class of diagrams involving 3 correlated particles. For the case k = 2, the simplest 3-particle diagram 
is the one where 3 particles leave a vertex, and bounce directly back on the subsequent advective step. This 
diagram gives \i = 1/27. Exactly solving the resulting equations for the 1-, 2- and 3-particle CCF's, we 
find that there are precisely 2 solutions (aside from the unstable solution at / = 1/2). Thus, it seems clear 
that the secondary equilibria generated by the 2-particle BBGKY equations are spurious, since they can be 
removed by such a simple perturbation. Unfortunately, including an arbitrary set of 3-particle diagrams, 
without performing the systematic 3-particle BBGKY approximation, tends to reduce the effectiveness of 
the approximation; thus, although the spurious equilibria are removed, the analysis described here does not 
give more accurate predictions for the actual equilibria than the 2-particle BBGKY analysis. To have a 
significantly improved approximation to the actual equilibria of the lattice gas, one would need to use a 
more complicated approximation scheme such as the complete 3-particlc BBGKY approximation. 

We conclude this section with a brief discussion of finite size effects. For any finite lattice, the complete 
equations of motion can have only a single equilibrium solution, corresponding to / = 1/2, since fluctuations 
can always drive a transition from one local equilibrium to another. Thus, if we have a lattice with I sites, 
the exact solution of the dynamical equations for all CCF's of 31 or fewer particles should only give a single 
solution. It is interesting to consider the effect that a finite lattice size would have on our discussion of 
the 2-particle BBGKY equations. The only way in which a finite lattice size would modify the equations is 
to change the coefficients 4>k(t) to correspond to random walks on the finite lattice. A particularly simple 
example of this is the degenerate case where we have a lattice with only a single vertex. In this case, the 
outgoing particles from a collision return immediately to the same vertex. Thus, we have 0fc(l) = 1 for all fc, 
and of course 4>k{t) = for all t > 1. This modification of the coefficients has no effect on the exact equations 
at a vertex, Eqs. ( p4[ ) and (p5f), so / = 1/2 is still a solution of the equilibrium equations. However, using 
the modified values for (f>, the 2-particle BBGKY equation (^) becomes 

C(/) = -4 + 42/ - 42/ 2 . 

This equation has two solutions, which give spurious equilibria analogous to those encountered previously 
on the infinite lattice. Thus, although the finite size effects remove the extra physical equilibria, which we 
only expect to exist in the thermodynamic limit, these effects leave the spurious solutions of the BBGKY- 
truncated equilibrium equations intact. An interesting question, which we will address in future work, is at 
precisely what lattice size the thermodynamic equilibria first appear in the 2-particle BBGKY approximation. 
An answer to this and related questions might shed light on the relationship between i-particle correlations 
and fluctuation scales. 



V. CONCLUSIONS 

We have described an NSDB lattice gas model for Schlogl's second chemical reaction. We derived a 
self-consistent set of equations for its exact homogeneous equilibria, solved these equations in the two- 
particle BBGKY approximation, and compared the results to numerical experiment. We found that this 
approximation describes the equilibria far more accurately than the Boltzmann approximation, but we also 
noted that it can give rise to spurious solutions to the equilibrium equations which can only be removed by 
including effects due to three-particle correlations. 

The possibility of the existence of spurious solutions of the two-particle BBGKY equations was raised by 
Bussemaker et. al.u The method they used to solve these equations was an iterative approximation method 
which was not well suited to recognizing the existence of multiple solutions. The use in this paper of a 
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diagrammatic formalism to describe the time development of the correlations made it possible to write the 
BBGKY-truncated equilibrium equations in a closed form which was amenable to numerical solution. It 
would be interesting to extend the diagrammatic analysis described here to higher-order truncations of the 
BBGKY hierarchy. 

The physically meaningful solutions of these BBGKY-truncatcd equilibrium equations provide an accurate 
description of the non-Gibbsian equilibrium of this lattice gas. The next step in this program of study will be 
to expand about this non-Gibbsian equilibria in Knudsen number, thereby generalizing the usual Chapman- 
Enskog analysis. In this way, the full reaction-diffusion equation, Eq. ([!]) will be derived, including the 
renormalized diffusion coefficient. This work is in progresstj. 
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TABLE I. Ensemble- averaged transition matrix 
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FIG. 1. The hexagonal lattice, with the checkerboard coloring and the enumeration of the three bits at each 
site. 
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FIG. 2. Evolution of the Schldgl model from random initial conditions yields domains of both low and 
high density, separated by sharp gradients whose width is governed by the diffusive term in the rate equation. 
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FIG. 3. Unique diagram contributing to 
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FIG. 4. Plot of <(/) versus / for k = 3. 
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Equilibrium Density 
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FIG. 5. Equilibrium density versus k. The black points with the error bars are from numerical experiment, the 
gray points without error bars are from the 2-particle BBGKY theory, and the line across the top is the Boltzmann 
value. 
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